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ABSTRACT 

We explore similarities and differences between several estimators of the cosmolog¬ 
ical bulk flow, B, from the observed radial peculiar velocities of galaxies. A distinction 
is made between two theoretical definitions of B as a dipole moment of the velocity field 
weighted by a radial window function. One definition involves the three dimensional 
(3D) peculiar velocity, while the other is based on its radial component alone. Different 
methods attempt at inferring B for either of these definitions which coincide only for 
a constant velocity field. We focus on the Wiener Filtering (WF, Hoffman et al. 2015) 
and the Constrained Minimum Variance (CMV, Feldman et al. 2010) methodologies. 
Both methodologies require a prior expressed in terms of the radial velocity correlation 
function. Hoffman et al. (2015) compute B in Top-Hat windows from a WF realization 
of the 3D peculiar velocity field. Feldman et al. (2010) infer B directly from the ob¬ 
served velocities for the second definition of B. The WF methodology could easily be 
adapted to the second definition, in which case it will be equivalent to the CMV with 
the exception of the imposed constraint. For a prior with vanishing correlations or very 
noisy data, CMV reproduces the standard Maximum Likelihood (ML, Kaiser 1988)) 
estimation for B of the entire sample independent of the radial weighting function. 
Therefore, this estimator is likely more susceptible to observational biases that could 
be present in measurements of distant galaxies. Finally, two additional estimators are 
proposed. 

Key words: Cosmology: large-scale structure of the Universe, dark matter 


1 INTRODUCTION 

On large scales where baryonic processes are dynamically 
unimportant, peculiar motions of galaxies are likely to be un¬ 
biased tracers of the peculiar velocity field of the dominant 
dark matter. This is in contrast to the spatial distribution of 
galaxies which is naturally biased with respect to the under¬ 
lying mass density field. Therefore, catalogs of radial pecu¬ 
liar motions should offer a unique window to any dynamical 
deviations from standard gravity which may arise in theo¬ 
ries for the observed cosmic acceleration (e.g. Hellwing et al. 
2014). One important caveat is that the extraction of cos¬ 
mological information from peculiar velocity catalogs gener¬ 
ally involves the biasing relation between galaxies and mass 
and also observational cuts imposed on the data, e.g. the 
Zone of Avoidance. An example is a direct calculation of 
the velocity correlation function where the desired signal is 
modulated by how galaxies are distributed in the particular 
velocity catalog (e.g. Davis & Peebles 1983). A large scale 
moment of the velocity field which has been the subject of 
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vivacious debate is the bulk flow . Although its estimation 
is affected by the spatial coverage of the data, the bulk flow 
is fairly insensitive to galaxy-mass biasing Li et al. (2012), 
Let v(r) be the three dimensional (3D) peculiar velocity as a 
function of the comoving distance coordinate r (in kms“^) 
and u{r) = v ■ r he the corresponding radial peculiar veloc¬ 
ity field. We consider two theoretical definitions of the bulk 
flow, B, 

B? = j <frg{r)v°‘{r) , (1) 

and 

B?i = 3 j A^rg{r)u{r)h°‘ {r) , (2) 

where h“ = are cosine angles between r and the Carte¬ 
sian axes dehned by unit vectors a:“ (a = 1, 2 & 3). The 
radial weighting function g{r) is usually introduced in the 
definition of B, and it satishes / d®gi(r) = 1. For a Top-Hat 
window of radius R centered on the observer g = 3/{A-kR^) 
for r > R and vanishes otherwise. For a Gaussian window 
g oc exp(—r^/2i?^). For a constant 3D peculiar velocity, 
V = Bo = const, the two definitions yield the same result. 
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i.e. Bq. However, they do not coincide in general (Nusser 
2014). 

Any moment of the peculiar velocity field can serve as 
a basis for a quantity that can be estimated from the data 
for the purpose of constraining cosmological models. The 
general framework for the various usages of the term bulk 
flow is the fact that the weighting function g is independent 
of direction. 

We aim at clarifying the connection between methods 
for inferring bulk flows from sparse and noisy velocity cata¬ 
logs. Different methods adopt one of the bulk flow definitions 
above (e.g. Feldman et al. 2010; Hoffman et al. 2015). The 
lack of equivalence between the two definitions even for a 
full 3D velocity field, guarantees different results for some of 
the methods. Further, although not always explicitly stated 
in the relevant papers, all methods rely on certain assump¬ 
tions on how the observed radial motions are related to the 
full 3D velocity field. The bulk flow is defined as integrals 
over space including regions uncovered by the sparse data. 
Therefore, even in the ideal case of no observational errors, 
the data on its own is insufficient to compute B. Supple¬ 
menting the missing information requires extrapolating the 
observed peculiar velocities to unobserved points in space. 
This is best done by resorting to the statistical nature of 
the 3D velocity field within the context of a cosmological 
model. A unwarranted criticism that is often raised is that 
the approach is circular in the sense that the inferred B is 
necessarily consistent with the assumed model. The target 
quantity (e.g. B) is certainly poorly constraint by very noisy 
data, in which case the prior information dominates. Nev¬ 
ertheless, current catalogs, e.g. the SFI-I—I- (Springob et al. 
2007) and Cosmicfiows-2 (hereafter CF2, Tully et al. 2013), 
are sufficiently accurate and dense and they constrain the 
bulk flow within ~ Mpc with very little dependence 

on the assumed prior (Nusser & Davis 2011). The methods 
discussed here will be addressed within the context of gaus- 
sian random fields. The assumption is justified since mo¬ 
tions on large scales obey linear theory for gravitational in¬ 
stability and hence, in the standard paradigm, the veloc¬ 
ity and density fields remain guassian. We further have to 
specify the statistical nature of the random error on the 
estimated peculiar velocity. Most common distance indica¬ 
tors are intrinsic scaling relations between galaxy observ¬ 
ables involving the distance modulus (log distance) rather 
than the actual distance. Thus a normal (gaussian) scat¬ 
ter in these relations propagates into a skewed distribu¬ 
tion of the measured distance (and peculiar velocity) er¬ 
ror (e.g. Lynden-Bell et al. 1988; Springob et al. 2014). The 
challenge of dealing with non-gaussian errors can be alle¬ 
viated by following the scheme of Nusser & Davis (1995, 
2011). Consider the (inverse) Tully-Fisher relation as a dis¬ 
tance indicator: g = sM — g + go + A 77 where s and go are 
constants, g is log the line-width, M and Ag is a random 
scatter with a gaussian distribution. For 2 <C 1, we have 
cz = r + u and M = m — 51og(r) = Mo — 51og(l — ujcz), 
where m is the apparent magnitude and Mo = m — 5 log(c 2 ) 
Typically, u/cz -C 1 even for nearby galaxies if u is mea¬ 
sured in Local Group frame. Nusser & Davis adopt A 77 = 
sMo + 2.17sulcz — g + go as the estimator for the veloc¬ 
ity where u is typically expressed in terms of a velocity 
model such as a bulk flow. This is equivalent to setting 
u = 0.46c2(r; — go — sMo)/s as the estimate of the pecu¬ 


liar velocity with a normally distributed random error. A 
closely related scheme has been advocated more recently by 
Watkins & Feldman (2015) for obtaining individual peculiar 
velocities. 

We also ignore here any systematic biases in the 
determination of the peculiar velocities. To mitigate 
spatial Malmquist biases (e.g. Lynden-Bell et al. 1988; 
Strauss & Willick 1995) we assume that galaxies are placed 
at their redshift coordinates rather than the observed dis¬ 
tance (Aaronson et al. 1982). To linear order this is equiv¬ 
alent to peculiar velocities expressed in terms of actual dis¬ 
tance and hence, for simplicity of notation, we shall assume 
that the velocities are give in terms of the actual distance. 

The outline of the paper is as follows. The notation and 
the statistical tools are given in §2. The “frequentist” ap¬ 
proach to inferring B is presented in §3. A new generalized 
ML estimation is presented here and the standard ML is 
shown to be recovered as a special case. In §4 we describe 
the relevant Baysian inference approach. §sec:MV focuses 
on minimum variance and constrained minimum variance. 
This section includes new estimator which incorporates con¬ 
straints in a probabilistic manner. A specific scheme for com¬ 
puting the relevant covariance matrices is given in §6. We 
end with a general discussion in §7. 

2 PRELIMINARIES 
2.1 Notation 

We are provided with a set of N data points di representing 
observations of the radial peculiar velocities of galaxies, 

di — Ui Gi , (^) 

where Ui is the underlying true signal and a is a ran¬ 
dom variable representing observational errors. In general 
(Ji=G^ + where ct* correspond to small scale veloc¬ 
ity dispersion and ad is proportional to the distance and it 
arises from the intrinsic scatter in the distance indicator (e.g 
the Tully-Fisher relation). The cosine angles, n“ = r ■ x°‘, 
between the radial direction and the Cartesian coordinates, 
satisfy the orthogonality condition, 

I , (4) 

where is the Kronecker delta function which equals to 
unity for a = /3 and vanishes otherwise. Let T7ij = (didj) 
and Sij = {uiUj) denote, respectively, the data-data and 
the underlying signal-signal correlation functions^. The an¬ 
gle brackets imply ensemble average over all possible real¬ 
izations of the quantity inside the brackets. Since the error, 
Gi, and the signal, m, are uncorrelated, we obtain 

T>ij = afS^J^ Sij . (5) 

The specification of the target quantity of interest, i.e. the 
bulk flow, is done via the cross correlation V, 

V? = (diB^^) = (uiR“) . (6) 

^ We shall use the terms “correlation functions” and “covariance 
matrices” interchangeably. 
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This is the only quantity which involves the definition of the 
bulk flow. For the definition (2) for example we get 

VT = S j rg{r)Si{r)h° (r) , (7) 

where Si{r) = (uiu{r)). It is also possible to obtain Vt 
for the first definition which involves the correlation with 
transverse component of v (Gorski 1988), but we do not 
present this here. 

We shall use boldface to denote vectors and matrices, 
e.g. d represents the explicit notation di for all i and "D is 
We will alternate between these two conventions, 
as demanded by the facilitation of mathematical manipula¬ 
tions. 


three component of the bulk flow and X 2 = d the N- 
dimensional data vector. The specification of the precise 
definition of the bulk flow is fixed via the covariance ma¬ 
trix Ei 2 = Esd = V = (uB). The matrix En = Ebb is 
a 3 X 3 matrix corresponding to {B°‘B^). Isotropy implies 
Egg = a%S^^. Further, E 22 = ^dd = "D. Therefore, accord¬ 
ing to (11), 

Ui = ^Ug2prB“ (13) 

a 

a 

and the minimization of the corresponding expression for A 
in (12) yields 


2.2 Probabilities 


The joint probability distribution function (PDF) P(B;d) is 
our main tool. Specifically, the conditional PDFs P(B|d) and 
P(<i|B) serve as the basis of the methods described here. We 
give a brief summary of the properties of conditional normal 
PDFs (e.g. Bertschinger 1987; Hoffman & Rlbak 1991). Let 
X be an n-dimensional gaussian random variable divided 
into two parts X = (Xi,X 2 ). The matrix E = (XX"'’) and 
its inverse, E“^ are decomposed as 


E = 


Ell 

S 21 


S 12 

S 22 


E"^ = 


rjiii 

j,21 




( 8 ) 


where Efej = {XkXj) and 


= Kb = (m) 

where 

= ^(<TB")"coVrpf. (15) 

ij 


3.1 A Special Case: The standard ML estimation 

As a special case we consider the velocity model v = Bq = 
const, implying Ui — nf ■ Bq . For either dehnition of B 
(1 & 2) and for any g{r), this model gives 


E“ = (Eii - Ei2E22iE?2) ' (9) 

= Ell + ^11 ^12 ^^22 — S12E11 Ei 2 ^ E12E11 

S'" = (E22-S?’2Sri'El2)“' 

= ^22^ +S22^Si2 ^Ell — El2E22^Si2j Ei2E22^ 

Ei= = -Eri^Ei2 (E22-S?’2Sn^Si2)'' 


Sij 

/ \ 2 \ ^ ^ q: /v q: 

= {UiUj} ^ as 2_^ni Uj 

(16) 

VT 

= {uiB°‘) = hfcr% 

(17) 



(18) 


— O'i Oij , 

(19) 


where the last equality is obtained by substituting (5) in the 
expression (13) for Bij. Therefore, (12) reduces to 


= (E-^ 


The joint PDF P{X) = P{Xi,Xi) ocexp(—Q/2) with 


A _ (di — J2a 

7. * 


( 20 ) 


Q = X?’Er/Xi + (X2 

-M)^r'(^2-M) 

(10) 

where 



/Lt = EJ 2 S//X 1 and 

^ = S 22 — E^2Si1^Si2 • 

(11) 


In this expression the residual between the data and the 
model for m is uncorrelated and is entirely due to observa¬ 
tional error and possible very small scale velocity dispersion. 
The solution for dA/dP“ = 0 is 


From this form of Q, the conditional PDF for X 2 given 
Xi is easily found using Bayes theorem, P(X 2 |Xi) = 
P(Xi,X 2 )/P(Xi) oc exp(—A/2) where 

A = (X2-M)^r'(^2-/i) • (12) 


P° 


N 

f3 i=l 



where 


( 21 ) 


3 INFERENCE BASED ON P(d|B): 
GENERALIZED ML ESTIMATION 

An estimate for B is obtained by means of Maximum Like¬ 
lihood (ML) estimation, i.e. by maximizing P(d|B), the 
likelihood for observing the data given the model. This is 
equivalent to finding B which renders a minimum in A by 
solving dA/dB°‘ = 0. We take Xi = B to represent the 




E 


( 22 ) 


The summation is over all particles independent of the form 
of g(r). Thus, this special velocity model reproduces the 
standard ML result (Kaiser 1988). In order to derive esti¬ 
mate for the B within a sphere of radius Rq, the actual 
model has to be modified into v(r) = Bo for r < Rq and 
■u = 0 otherwise. 
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4 INFERENCE BASED ON P(B|rf): WIENER 
FILTERING 


4.2 Equivalence with B from a reconstruction of 
the full 3D peculiar velocity field 


We now consider the conditional PDF P(B|d). Therefore, 
we take Xi = d, X 2 = B, En = 2? and, as before, E 12 = V. 
The corresponding expression for A is A = (B — — 

y.) and its minimum is rendered at 

= (23) 

where the last equality is obtained by substituting the rel¬ 
evant quantities in (13). The covariance matrix ^ does not 
appear in B^i^, but we give it here for later use 

- Y, Vfpf . (24) 


4.1 The relation between B^|^ and B^|^ 

Since P(B|</)P(d) = P(<i|B)P(B), the expressions (23) and 
(14) must coincide in the limit of P(B) = const, i.e. for 
(jb —> 00 . To see how this happens we resort to the matrix 
notation and write A°‘^ in (15) as 

A = EggEsd (^dd — ^Bd^BB^Bd^ ^B<i^Ss(25) 

E bb 1 

~ ^BB ? 


where the equalities are derived using (9) with 1 —>■ P and 
2 —>■ d. Further the sum over i,j in the solution (14) trans¬ 
forms into 


-‘BB^Bd 


■-Bd^BB^Bd 


^ ^dd ^Bd I Sbb — Er^E.,., E 


Bd^dd ^Bd 




= ^dd ^Bd 
= 

Hence, in matrix notation (14) is 

With a little more matrix manipulation we obtain 


'd\B ~ ^ ^B\d 


B^ib = 


Sbb — 






(26) 


(27) 

(28) 


where the last equality is obtained from the first two lines 
in (9). The contribution of (E^®) is negligible compared 
to Ebb in the limit of ctb —>■ C) 0 . This is demonstrated by 
proving. 


Hoffman et al. (2015) use the WF methodology 
(Zaroubi et al. 1995) to derive the full 3D peculiar ve¬ 
locity field in space, v{r), from the measured radial peculiar 
velocities in the CF2 catalog (Tully et al. 2013). They then 
compute the bulk flow as defined in (1) for a Top-Hat 
window, i.e. the mean v within a sphere centered on the 
observer. They could have also computed the bulk flow 
according to definition (1). We will show now that B derived 
from this procedure is (unsurprisingly) entirely equivalent to 
B^ij. The full field is obtained as the one which maximizes 
the probability P{v{r)\d). Substituting Xi = d and X 2 = v 
in (11) and minimizing A = —21nP -|- const in (12), we 
obtain 

u“(r) = ^p-4(niu“(r))d, . (31) 

Substituting this expression into either of the bulk flow def¬ 
initions (1) & (2), leads to the same expression as B^|^ in 
(23). 

5 MINIMUM VARIANCE ESTIMATION 

So far we have considered estimates which maximize normal 
conditional PDF. The minimum variance (MV) approach 
described below yields identical results, but is not restricted 
to normal PDFs (cf. Zaroubi et al. 1995, for a detailed ac¬ 
count). A MV estimate of B is sought as a linear combina¬ 
tion of the data, 

3 N 

B^^v = w’c* = E ■ (32) 

q: = 1 i=l 

The weights Wi are found by minimizing the variance, Amv, 
of the residual between Bj^,^ and all of its possible realiza¬ 
tions, 

A^^ = ((B-md)^) (33) 

= (B^) — 2Vw + vADw . 

The minimum in Amv is obtained at 

w'mv = ■ (34) 

Implementing this result into (32) yields = 

j Therefore, referring to (23), we find B^jy = 

Bfiid- 


5.1 Constrained Minimum Variance 


(B'^ (e^®) ' B) < (B^EbbB) = 3a| . (29) 

We write di = yi-\-/S.i where yi is the mean value of di given 
B (cf. eq. 13) and Ai is a random variable uncorrelated with 
P“. Hence, in the limit of very large ctb, 

= (30) 

a 

Substituting (E®®) ^ = Ebb — ERdS))'j^Esd correspond¬ 
ing to , the l.h.s of the inequality 

becomes Sa% - eI .;3 Eij (B“P''} « 0, for the 

approximate T> given in (30). 


One may impose constraints to any of the above proce¬ 
dures, as done by Feldman et al. (2010) for the MV esti¬ 
mator. Their constraint is that if the 3D velocity field is 
V — Bo = const, then (32) must reproduce Bq in the absence 
of observational errors. Substituting di — m = Ec Bq Af in 
(32), the constraint implies 

= (35) 

i,l3 

must hold for any Bo = const, yielding 

N 

E ■ (36) 

i=l 
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This constraint is incorporated by modifying (33) into 

Acmv = ((B“ - E + E . 

i 13 \ i J 

(37) 

to be minimized with respect to the Lagrange multipliers, 
A“^ as well as wf. The solution for the minimum point is 
(c.f. eqs. 8-10 in Agarwal et al. 2012), 




(38) 


3 

= E 

7^1 - 

where 

""' = )Eo5‘’ 




(39) 


Feldman et al. (2010) adopt the following scheme to 
approximate (in §6 where we outline an alternative 
scheme). They write 


N' 


{diB°') = E w'fidiUji) 
j'=i 

where 


(40) 


(41) 


are the weights of an isotropic survey of N' exact radial 
velocities «;/ measured at random positions r[, having the 
radial distribution fixed by g{r), and 


N' 

j\r/ / ^ “ * 


(42) 


In the continuum limit of A^' oo, = 

dtlBdrg(r)h°‘(r)n^(r)/ J dO,r^drg{r) = l/3(5((^. Hence, 
w'ii = 3n“ and (40) becomes 

{diB°‘)=3 J dQr'‘^dr'g{r')n°‘{r){uiu{r)), (43) 

which clearly coincides with (7). 


where 




E 


Therefore, 


w = w 


ML 5 


(47) 


(48) 


i.e. the weights reduce exactly to the ML weights. In con¬ 
trast, Bg|j approaches zero for uncorrelated velocities, as 
expected. 


5.3 Generalization of the constraint 

The constraint imposed above yields weights that recover 
a constant velocity field in the case of perfect data. The 
data points are not uniformly distributed and one may argue 
that this constraint is actually inappropriate since there will 
always be leakage from other modes of the velocity held. 

If desired, other constraints which may be more real¬ 
istic could be imposed. In particular, the radial velocity 
held could be decomposed into a functional basis of prod¬ 
ucts of spherical harmonics and spherical Bessel functions 
(Regos & Szalay 1989; Scharf & Lahav 1993; Fisher et al. 
1995) . If the velocity held is fully specihed by a dipole term 
then we could write u{r) = CL°'‘>p{r)h{r)°' where a“ are 

the expansion coefhcients and ip represent the radial func¬ 
tional basis, e.g. a Bessel function with wavenumber k (see 
§6 for details). The constraint is that for the ideal case of this 
form of the velocity held, the bulk how is correctly recov¬ 
ered. Adopting the second B dehnition in (2), the constraint 
becomes 

3 f d^rff(r) a^ip{r)n^{r)n°‘{r) 

= E E (49) 

z P 

Performing the angular integration and remembering that 
the constraint must hold for any we get 

E (50) 

i 

where d' = 47r J drr^g{r)ip{r). Substituting this equality as 
our constraint in Acmv (instead of 36), the corresponding 
weights are given by (39) but with the modihcation, hf —>■ 
iPifiT and (5“'^ 


5.2 No velocity correlation &z standard ML 

We discuss now the case where Dij = (crj +aQi)Sij . Accord¬ 
ing to (40), {diB°‘) = 0 since overlap between the data and 
the ideal sample is zero. Therefore, 


■\joiB 1 V ^ —2^ol B 

M = — > cr. n,- m 

(44) 

1 


= 


= -(M-i)“^ 

(45) 

07 

<5 

CO 



(46) 


5.4 Alternative method for imposing constraints 

Constraints can be imposed by modifying the conditional 
PDF P(B|d) oc exp(—A/2). This is done adding either a 
quadratic or linear term in B to A. Here we present only 
the quadratic modihcation since it leads to a particularly 
elegant result. We write 

A„,od= (B-B^„)^r'(B-B^„)+BT(4A)-iB. (51) 

as minus twice the log of the modihed PDF, where B^|^ 
and ^ are given in (23) and (24), respectively. The 3x3 
matrix, A, is derived by demanding that B which renders 
a minimum in Amod also satished the desired constraint. 
For the constraint of Feldman et al. (2010) of B = Bq for 
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di = f • Bo we find (A D^■^Vfn^ and the 

solution is given by, 


WB = B 

** -‘-'mod ^B\d ’ 


(52) 


where Since = Ei«'Mvi'^i- the 


result clearly reproduces Bo for di = Ec hf Bo • 


6 COMPUTING THE COVARIANCE 
MATRICES 

In a cosmological model with gaussian initial conditions, all 
statistical properties of the linear density contrast field, 5 are 
fully specified by the linear density power spectrum, P{k), 
defined via = (27r)®5 {k — k')P{k) where Sj. is 

the Fourier transform of S{r) and k is the wavenumber. As¬ 
suming potential flow, the velocity statistics can readily be 
related to the density field through the velocity-density re¬ 
lation S = —/(r2in)V ■ V (Peebles 1980), where f{Qrn) is the 
linear growth factor and flm is the mass density parameter. 
All relevant correlation functions can be expressed in terms 
of definite integrals involving these power spectra. The inte¬ 
grals must be evaluated numerically with lower and upper 
bounds on k corresponding to small scales out to horizon 
scales. The integrals are typically cumbersome with highly 
oscillatory functions (spherical Bessel) appearing in their in¬ 
tegrands. Instead of direct numerical integration, we propose 
evaluating the integrals using the formalism developed by 
Fisher et al. (1995) and this section heavily relies on that 
paper. Radial peculiar velocities are conveniently expressed 
in terms of spherical Bessel functions, ji , and spherical Har¬ 
monics, Yim- To do that we imagine a very large spherical 
volume or radius Bmax which entirely encompasses the data 
and is also sufficiently large that the effects boundary con¬ 
ditions are small. Expressing the density contrast as 

^max ^max 

^w=EE E Clndlmnjl{knAYlm{^') (^^) 

Ti—1 l—O m— — l 

where the wavenumbers are fixed by the boundary condi¬ 
tions at Bmax (Fisher et al. 1995). We advocate the bound¬ 
ary conditions which yield gravity potential decaying as r“* 
for r > Bmax. Since jiYim are eigenfunctions of the Lapla- 
cian, these boundary conditions are equivalent to finding kn 
which satsify d\nji{k„r)/d\nr\R^^^ = —{I + 1). The num¬ 
bers Cin are fixed entirely by the boundary conditions and 
are given in Table Al of Fisher et al. (1995). The density 
coefficients are 

Simn = f d^rS(r)ji{knr)Yi^{r) . (54) 

Using the linear velocity-density relation one obtains 
(Regos & Szalay 1989; Scharf & Lahav 1993; Fisher et al. 
1995), 

u{r) = fYl Yimir) . (55) 

, 

l,m,n 

where ji{z) — Aji[z )/For velocities measured 
with respect to the Local Group frame the I — 1 term should 
be modified by subtracting 1/3 from j[{knr). Using 

((5qmini(5r2m2n2> = P > (56) 


this representation for u{r) yields 

5(n,r,) (57) 


h.2 

, 

l,n 

= E -^-^^fi(^r,ri)j',{k„rj)Pi{ri ■ rj) 

l,n 

where Pi is Legendre polynomial of order I and we have 
used 47r/(2Z -|- 1) E^ Yim{ri)Y*^{rj) = Bi(fi ■ Tj). From S, 
the cross correlation can easily be computed P. Using the 
expression (7) which is appropriate for for definition (2) we 
get^ 




(58) 


3fhT 


P{kn)jl{knri) 

kl 


ArA g{r)j[{knr) 


7 DISCUSSION 

An aim of this contribution is to show that different ap¬ 
proaches to the determination of B are tightly related. All 
methods yield something which mimics the behavior of the 
theoretical bulk flow. The non-trivial challenge is to con¬ 
trast any estimate with predictions of cosmological model. 
Systematic errors related to observables of galaxy properties, 
contaminate various methods in different ways. More distant 
galaxies are more prone to these systematics, therefore, con¬ 
sistency checks must be performed by applying the methods 
subsamples of the data. Furthermore, spatial coverage of the 
data depends on galaxy-type and unavoidable cuts on the 
observables. Thus, methods should be tested on mock galaxy 
catalogs that mimic the observations in as much as possible, 
including the observed large scale structure. The constrained 
simulations (Sorce et al. 2013) designed to match the low 
redshift large scale structure, can potentially be highly ben¬ 
eficial if combined with galaxy formation models. 

We have also emphasized the differences between the 
methods. Hoffman et al. (2015) provide an estimate for the 
bulk flow as defined in (1), while Feldman et al. (2010) aim 
at B given by (2). These two definitions of B do not coin¬ 
cide (Nusser 2014). and should not produce identical results. 
Nonetheless, both methods can easily be adapted to any of 
the two definitions. 

In §3 we derive a new B estimator which naturally 
arises from the likelihood B(rf|B). This is a generalization of 
the standard ML of Kaiser (1988), but incorporates veloc¬ 
ity modes beyond the (constant) bulk flow of the entire data 
catalog. Essentially, this estimator is equivalent to marginal¬ 
izing P{v^^^^^\d)/PCB) over all modes of a general velocity 
model, i excluding a mode corresponding a global bulk 

flow . Another new estimator is given in §5.4. This estimator 
incorporates the Feldman et al. (2010) constraint in a new 
and simple way. But we emphasize that in a Baysian ap¬ 
proach, the prior as expressed in S and P correlation func¬ 
tions. already contains all missing information. Additional 

^ The following integrals have been used J dUh^Yio = 
- /dUh-Un = , /dUh-yi_i = , 

/dOn^Ull = and /dUh^U-i = 
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constraints (not based on extra information) are either sta¬ 
tistically redundant or incompatible with the model prior. 
For example, requiring the estimator to produce Bq = const 
for the input data di — Vi ■ Bo (no observational errors), ig¬ 
nores variations in the velocity field on the scales that are 
not probed by the data. 

We have refrained from making any quantitative assess¬ 
ment of the difference between methods and the bulk flow 
definitions. We hope to do that in the near future as well 
as comparison applying all methods to observational datat, 
but it is not the point of the paper. 
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